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ABSTRACT 

Our research objective in this paper is to reconstruct an initial linear density field, 
which follows the multivariate Gaussian distribution with variances given by the linear 
power spectrum of the current CDM model and evolves through gravitational instabil- 
ity to the present-day density field in the local Universe. For this purpose, we develop 
a Hamiltonian Markov Chain Monte Carlo (HMC) method to obtain the linear density 
field from a posterior probability function that consists of two components: a prior of 
a Gaussian density field with a given linear spectrum, and a likelihood term that is 
given by the current density field. The present-day density field can be reconstructed 
from galaxy groups using the method developed in Wang et al. (2009). Using a realistic 
mock SDSS DR7, obtained by populating dark matter haloes in the Millennium sim- 
ulation with galaxies, we show that our method can effectively and accurately recover 
both the amplitudes and phases of the initial, linear density field down to a scale of 
about 2.8/i" 1 Mpc. To examine the performance of our method in the highly non-linear 
regime, we use ./V-body simulations to evolve these reconstructed initial conditions to 
the present day. The resimulated density field thus obtained accurately matches the 
original density field of the Millennium simulation in the density range 0.3 ^ p/p ^ 20 
without any significant bias. Especially, the Fourier phases of the resimulated density 
fields are tightly correlated with those of the original simulation down to a scale cor- 
responding to a wavenumber of ~ 1 ftMpc -1 , much smaller than the translinear scale, 
which corresponds to a wavenumber of ~ 0.15 /iMpc -1 . 

Key words: dark matter - large-scale structure of the universe - galaxies: haloes - 
methods: statistical 



1 INTRODUCTION 

In the current CDM cosmogony, a key concept in the build- 
up of structure is the formation of dark matter halos. 
These are quasi-equilibrium systems of dark matter, formed 
through non-linear gravitational collapse. In a hierarchical 
scenario like CDM, most of the mass at any given time is 
bound within dark halos; galaxies and other luminous ob- 
jects are assumed to form by the cooling and condensation 
of baryonic gas within these halos (see Mo, van den Bosch & 
White 2010). With current iV-body simulations, the proper- 
ties of the CDM halo population, such as the mass function, 
the spatial clustering, the formation history, and the internal 
structure are well understood. However, how galaxies form 
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in dark matter halos in the cosmic density field remains 
an unsolved problem. A long-standing problem in current 
galaxy formation theory is to explain the low efficiency with 
which baryonic gas is converted into stars: the observed mass 
in stars at the present time is less than 10% of the total bary- 
onic mass in the universe (Bell et al. 2003). Including cold 
gas associated with galaxies only increases this to ~ 12%. 
This low efficiency of star formation and gas assembly into 
galaxies is not a natural consequence of hierarchical forma- 
tion, in which the gas is expected to cool rapidly at high 
redshift in low-mass dark matter halos. A number of physi- 
cal processes have been proposed to suppress gas cooling and 
the star formation efficiency. These include photoionization 
heating by the UV background (e.g. Efstathiou 1992; Thoul 
& Weinberg 1996; Somerville 2002; Gnedin 2000; Hoeft et 
al. 2006), feedback from supernova explosions (e.g. White & 
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Rees 1978; Dekel & Silk 1986) and from AGN (e.g. Tabor 
& Binney 1993; Ciotti & Ostriker 1997; Hopkins et al. 2006 
and references therein; Wang et al. 2012a), and pre-heating 
by star formation/ AGN (e.g. Mo & Mao 2002; Pfroemmer et 
al. 2012), and by pre-virialization (Mo et al. 2005). Unfortu- 
nately, our understanding of all these processes is still poor, 
making it difficult to test the predictions of these scenarios 
with observations. 

A key step in understanding the galaxy formation pro- 
cesses throughout the cosmic density field is to study the 
distributions and properties of galaxies and the intergalac- 
tic medium (IGM), and their interactions with each other 
and with dark matter. In the local universe, detailed ob- 
servations of the galaxy population are now available from 
large redshift surveys, such as the Sloan Digital Sky Survey 
(SDSS; e.g. York et al. 2000). This not only allows us to 
study the large-scale structure in the local universe, but can 
also be used to derive a large number of physical quantities 
characterizing the intrinsic properties of individual galaxies, 
such as luminosity, stellar mass, color, morphology, size, star 
formation rate, nuclear activity. There have also been ob- 
servational programs dedicated to the various aspects of the 
IGM. Extensive X-ray observations have been conducted to 
study the hot gas associated with clusters and rich groups of 
galaxies but one expects the total mass associated with such 
systems to be small. With the advent of accurate measure- 
ments of the cosmic microwave background from observa- 
tions such as the South Pole Telescope (SPT), the Atacama 
Cosmology Telescope (ACT), and the PLANCK Satellite, 
one can also probe the hot, diffuse gas outside clusters and 
groups through the Sunyaev-Zel'dovich (SZ) effect. However, 
as shown in Mo & White (2002), at the present time about 
70% of all the mass is in virialized halos with virial tem- 
peratures below 10 6 K, too cold to be studied with X-ray 
data and/or the SZ effect. A promising way to study the 
diffuse IGM at such low temperature is through quasar ab- 
sorption lines. With the installation of the Cosmic Origins 
Spectrograph (COS) on the HST the sample of UV absorp- 
tion systems at low redshift is expected to increase by an 
order-of-magnitude or more, allowing a much more detailed 
examination of the warm component of the local IGM. 

These observational programs together provide an un- 
precedented data base to study how galaxies form and evolve 
in the cosmic density field. However, in order to make full 
use of the potential of the observational data to test models, 
one has to develop an optimal strategy. Conventionally, one 
starts with a cosmological model, e.g. the current ACDM 
model, uses computer simulations to follow the evolution 
of the cosmic density field, and compares simulation results 
with observational data in a statistical way. However, such a 
comparison can be made directly only under the assumption 
that the observational sample and simulation are fair rep- 
resentations of the universe, so that cosmic variance is not 
an issue. Unfortunately, this assumption is almost always 
violated in reality. Simulations are limited by the dynami- 
cal ranges they can cover. In order to resolve processes on 
the scale of galaxies, the simulation volume has often to be 
much smaller than a fair sample of the large scale structure. 
Observationally, finite sample volumes also lead to a biased 
representation of the statistical properties of the large-scale 
structure and galaxy population in the universe. 

It is, therefore, imperative to have as much theoreti- 



cal and empirical input as possible both to design an opti- 
mal observational strategy and to help interpret the limited 
amount of observational data in an unbiased way. The un- 
certainties are minimized if comparisons between observa- 
tion and model prediction are made for systems that have 
both the same environments and the same formation his- 
tories. Ideally, if we can accurately reconstruct the initial 
conditions for the formation of the structures in which the 
observed galaxy population resides and from which the ac- 
tual gas emissions and absorptions are produced, then we 
will be able to compare observation and simulation (i.e., 
data and theory) in a way that is free of cosmic variance, 
thereby greatly enhancing the constraining power of the ob- 
servational data. 

The goal of this paper is to develop a method that can 
be used to reconstruct the initial (linear) density field that 
forms the density field in the local universe. In this first 
paper in a series, we describe our reconstruction method and 
test its performance with realistic mock galaxy catalogs. The 
structure of the paper is arranged as follows. In Section [2] 
we describe our reconstruction method. In Section^ we test 
our method using a simulated density field. In Section!?] we 
present our mock catalog that is used to test our method and 
our results of the linear density field reconstructed from it. In 
Section [5] we use N-body simulations to follow the structure 
formation seeded by the re-constructed initial density field, 
and compare the final density field with the original density 
field used in the mock catalog. Finally Section [6] contains a 
summary of the main results and some further discussions. 

In order to avoid confusion, we here list the various 
matter density fields used in the text: (i) The final den- 
sity field, pf(x): the true present-day density field, either in 
the real Universe or in a cosmological iV-body simulation, 
(ii) The reconstructed density field, p rc (x): the present-day 
density field re-constructed from the mock catalog, (iii) The 
reconstructed initial (linear) density field: the initial linear 
density re-constructed from a present-day density field, (iv) 
The re-simulated density field, p rs (x): the present-day den- 
sity field obtained from numerical iV-body simulations using 
the re-constructed initial density field as initial conditions, 
(v) The modeled density field, p mo d(x): a model prediction 
for the present-day density field obtained from the initial 
density field using the modified Zel'dovich approximation 
introduced by Tassev & Zaldarriaga (2012). This model den- 
sity field is used in the reconstruction method to link the 
initial and final density fields. 



2 METHOD 

Our reconstruction consists of the following steps: (i) Use 
galaxy group^] selected from redshift surveys, such as the 
SDSS, to represent dark matter halos; (ii) Use dark mat- 
ter halos above a certain mass to reconstruct the cosmic 
density field at the present day; (iii) Reconstruct the initial 
density field that best matches the final density field under 
the constraint of current cosmology and a linear perturba- 
tion spectrum. The galaxy group finder used is described 

1 We use galaxy group to refer a galaxy system (a cluster or a 
group) without regard to its richness. 
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and tested in detail in Yang et al. (2005, 2007). The method 
for reconstructing the density field starting from a sample 
of dark matter haloes (i.e., galaxy groups) above a given 
mass threshold is described and tested in Wang et al. (2009, 
2012). In what follows we describe how we reconstruct the 
initial linear density field from a given present-day density 
field. 



2.1 Objectives and the Posterior Probability 
Distribution 

Our goal is to obtain the linear density field that can repro- 
duce a given present-day density field. We work in Fourier 
space, so that the initial density field is specified by <5(k), the 
Fourier modes of the initial density field. Two constraints are 
used in the reconstruction. First, according to the standard 
cosmology, we assume the linear density field to be Gaussian, 
so that the Fourier modes obey the following probability dis- 
tribution: 



k j=0 



Trfltot*)] 1 /: 



.exp{-[5 J (k)] 2 /Piin(fc)} 



(1) 

where Pn n (k) is the (analytical) linear power spectrum, and 
the subscripts j = 0, 1 denote the real and imaginary parts, 
respectively. Because <5(k) is the Fourier transform of a real 
field, <5(k) = 8*(— k) so that only the Fourier modes in the 
upper half-space (i.e. with k z > 0) are needed. Second, the 
density field, p mo d(x), evolved from this linear density field 
according to a chosen model of structure formation, should 
best match the present-day (final) density field, pt(x). In 
other words, we seek the appropriate <5(k) to minimize a 
'cost parameter' which we define as 



E 



[Pmod(x) - pf(x)] 2 q;(x) 
2a t 2 (x) 



(2) 



where <Tf(x) is the statistical uncertainties in the final, 
evolved density field, Pf(x), while w(x) is a weight func- 
tion used to account for the survey geometry. As mentioned 
above, the final density field can be reconstructed from, 
for example, a galaxy redshift survey, with the method de- 
veloped in Wang et al. (2009) and we will come back to 
this in Section \4. 2 1 The uncertainties 0"f(x) are found to be 
roughly proportional to pf(x) (see Section [4.2|l . and so we 
set <Jf(x) = /xpf(x), with \x a constant parameter. In order 
to obtain the model prediction for the final density field, we 
need a model to link p mo d(x) with <5(k). This model should 
not only be accurate, but also be efficient so that the com- 
putation can be achieved in a reasonable amount of time, 
as to be discussed in Section \2. 31 In practice, all these fields 
are to be sampled in a periodic box of length L on a side, 
divided into iV c grids in each dimension, so that the number 
of Fourier modes to be dealt with is finite. 

Because of the statistical uncertainties in pf (x) and the 
finite survey volume, the solution for <5(k) under the two 
constraints described above is not unique. What we can 
seek is therefore a sample of solutions described by the pos- 
terior probability distribution of <5(k). Assuming that the 
likelihood of p mo d(x) given Pf(x) is exp(— x 2 ), the posterior 
probability distribution for <5(k) given pt(x) can be written 



Q(<5 J (k)|p f (x)) = e- s - [ "'«° d(x) -« (x)l2 " (x)/2tT f (x) 
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-[5 3 (k)] 2 /P Iln (fc) 
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In general sampling the posterior distribution in a high 
(of order N^) dimensional space is computationally chal- 
lenging. Here we propose to use the Hamiltonian Markov 
Chain Monte Carlo technique (HMC, Duane et al. 1987; 
Neal 1996), which has proven to be effective for explor- 
ing large, multi-dimensional posterior spaces (e.g. Hanson 
2001). Different from the conventional Markov Chain Monte 
Carlo, the HMC method introduces a persistent motion of 
the Markov Chain when exploring the parameter space so 
that the random walk is greatly suppressed and the effi- 
ciency much improved (Duane et al. 1987). This method 
has already been widely applied in astrophysics and cos- 
mology (see Hajian 2007; Taylor, Ashdown & Hobson 2008; 
Jasche & Kitaura 2010; Kitaura et al. 2012; Jasche & Wan- 
delt 2012, hereafter JW12). Once convergence is achieved, 
one can draw posterior samples to study the statistical prop- 
erties of the target distribution. In general, this kind of sta- 
tistical analysis requires a careful design of the likelihood 
function to take into account the detailed statistical uncer- 
tainties in the observational data. However for our purpose 
the requirement is less stringent. Indeed, as discussed above, 
our goal is to get the initial linear density field which obeys 
the probability distribution Eq. ([T} and will evolve under 
gravitational instability into a final state that is similar to 
the given present-day cosmic density. As we will demon- 
strate below using realistic mock catalogs constructed from 
cosmological TV-body simulations, our HMC method based 
on the likelihood function defined above is sufficient for our 
research objective. 



2.2 The Hamiltonian Monte Carlo Method 

In this subsection, we describe the Hamiltonian method in 
more detail (see also Hanson 2001; Taylor et al. 2008; JW12). 
The method is itself based on an analogy to solving a phys- 
ical system in Hamiltonian dynamics. As a first step, we 
define the potential of the system to be the negative of the 
logarithm of the target probability distribution, 

tf(*,(k)) = -ln[Q(«i(k)|pf(x))] 

half half 1 
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For each 5j(k), a momentum variable, f>j(k), and a mass 
variable, rrij (k) , are introduced. The Hamiltonian of the fic- 
titious system can then be written as 



half 1 2/. \ 



k 3=0 



2m f (k) 



(5) 



The statistical properties of the system is given by the parti- 
tion function, exp(— H) , which can be separated into a Gaus- 
sian distribution in momenta pj (k) multiplied by the target 
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distribution, 

half 1 pj(k) 

exp(-ff) = Qfc(k)|p f (x)] HHe^J^ . (6) 

k 3=0 

Thus, the target probability distribution can be obtained by 
first sampling this partition function and then marginalizing 
over momenta (i.e setting all the momenta to be zero). 

In order to sample from the partition function, we first 
pick a set of momenta pj(k) randomly from the multi- 
dimensional un-correlated Gaussian distribution with vari- 
ances rrij(\i). We describe how to pick the mass variables 
in Section \2. 41 We then evolve the system from the starting 
point [<5j(k),Pj(k)] in the phase space to some pseudo time 
T according to the Hamilton equations, 



dMk) 



dH 



9p 3 -(k) mj(k) 



Pj( 



d Pj (k) 
dt 



dH 



28, (k) 



88, (k) P lin (fc) 



(7) 



(8) 



where Fj(k) = dx 2 /d8j (k) is the likelihood term of the 
Hamiltonian force to be discussed below. The integrated tra- 
jectory finally reaches a point [<^(k),p^(k)] in phase space 
and we accept this state with a probability 



{l,« 



- IHty (k) ,j>$ (k))-H(^ (k) , Pj (k))] ' 



(9) 



The procedure is repeated by randomly picking a new set of 
momenta. 

Since the Hamiltonian of a physical system is conserved, 
the acceptance rate should in principle be unity, which is one 
of the main advantages of working with the partition func- 
tion, exp(— H), instead of the target distribution function 
itself. However, rejection may occur because of numerical 
errors. In order to optimize the acceptance rate, it is com- 
mon practice to integrate the "equations of motion" using 
the leapfrog technique, 



(10) 



8, (k, t + t) = *j(k, t) + — 7rrPi(k, * + r/2) ; (11) 
rrij (k) 



w (M + T )= a CM + T/2)-I^ 



(12) 



where r represents the time increment for the leapfrog step. 
The leapfrog technique uses half-time steps in the first and 
third equations so that the scheme is accurate to second or- 
der in r. In addition, this technique has the advantage of 
exact reversibility to ensure that the chain satisfies detailed 
balance (Duane et al. 1987). The equations are integrated 
for n steps so that nr = T. The value of T must be ran- 
domized to avoid resonance trajectories. We thus randomly 
pick n and r from two uniform distributions in the range 
of [l,n max ] and [0,r max ], respectively. We will discuss our 
choices of n max and r max below. The n leapfrog steps are 
referred to as one chain step. 



2.3 Hamiltonian Force and Structure Formation 
Model 

As shown in Eq. (JSJ), the Hamiltonian force consists of two 
components, the prior term, 2<5j(k)/Pii n (fc), and the likeli- 
hood term, Fj(k). The later can be re- written as: 

[Pmod(x) - Pf(x)]w(x) <9p mo d(x) 



*i(k) = £■ 



a 2 (x) dS 3 (k) 

9pmod(x) 



}2 Pd(k2)e * 2 " 9^k) Z) P m °d(ki)e* i x 



lV 3 X^n*(h- ^Pmod(kl) 

iVc^Pd(kx) g5j(M) 



(13) 



For the sake of simplicity we have introduced in this equation 
a new quantity Pd(x), which is directly related to p mo d(x), 
as defined in the second equation. 

To derive the solution of the Hamiltonian force, we 
need a model of structure formation to connect p mo d(k) 
and <5(k). In this work, we adopt the model developed in 
Tassev & Zaldarriaga (2012, hereafter TZ12). According 
to TZ12, a fully evolved density field, pf, can be written 
in terms of the modeled density field, p mo d, as Pf(k) = 
Pmod(k) + p mc (k). Here, p mc (k), a random mode-coupling 
residual, is generally small on mildly non-linear scales and 
can be neglected. The modeled density can be obtained via 
Pmod(k) = #,s(fe)pMz(k), where pmz is the density field 
based on the Modified Zel'dovich approximation (MZA) de- 
veloped by TZ12 and Rs{k) is a density transfer function. 
The transfer function is obtained via comparing the predic- 
tion of the Zel'dovich approximation and the A r -body sim- 
ulation, and so are the other transfer function shown be- 
low (see TZ12 for details). As demonstrated in TZ12, the 
MZA can predict the particle positions and velocities even 
slightly better than second-order Lagrangian perturbation 
theory. Following TZ12, we use the MZA to derive the dis- 
placement field, s(q), 

s(q) = ^s(k 2 )e* 2 ' q = ^7? z (fc 2 )s z (k 2 )e* 2q 
k 2 k 2 
ik 2 



= ^i?,(fe) T #<5(k 2 )e i 

k 2 



k 2 q 



(14) 



where q is a Lagrangian coordinate, s z (k 2 ) is the Zel'dovich 
displacement field, and J? z (fe 2 ) = exp(— 0.085(fc2/fcNL) 2 ) is 
the transfer function for the Zel'dovich displacement field 
with /cnl the non-linear scale at redshift zero. We move 
particles, which are initially located on uniform grids of 
positions q, to x p (q) = q + s(q) to sample the den- 
sity field. We then utilize a clouds-in-cells (CIC) assign- 
ment (Hockney & Eastwood 1981) to construct the MZA 
density field on grids from the particle population. We 
Fourier transform the MZA density, and multiply it with 
the density transfer function R$(k) = exp(0.58d), where 
d = <5 2 (fc/2),0 and a Gaussian kernel wo{R B k) character- 



2 Both k NL and <5 2 can be read off from Eq. (3.9) in TZ12. The 
detailed form of the density transfer function for MZA was com- 
municated to us by Svetlin Tassev. 
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ized by a smoothing scale R s . We then deconvolve the CIC 
kernel by dividing the resulting density field in Fourier space 
by the Fourier transform of the CIC kernel, u>cic(k) = 
smc(k x L/2N c )smc(k y L/2N c )smc(k z L/2N c ), where k x , k y 
and fc z are the x, y and z components of the wavevector 
k. Finally we obtain the modeled density field as, 

Rs(ki)w G (Rski 



Pmod 



(ki) = 



^E< 



-ikiXp(q) 



(15) 



Inserting Eqs. (|14|l and (|15[1 together with x p (q) = q + 
s(q) into Eq. (|13|) then yields the likelihood term of the 
Hamiltonian force: 



^(k) = 



k! x p (q) gHki ■Xp(q)) 



x E e 

q 



E 



R z (k 2 ) dSfa) 
k\ dSjQt) 



(*»•)£< 



ik 2 q 



x^(-jki)pS(ki)i? i (fci) U ;G(i?Bfci)e- lkl ' x >' (q) . 

ki 

(16) 

Note that dx p /d5j(k) — ds/dSj(k) is used in the derivation. 
For convenience we introduce a density-vector field, 

*(q) = N? ^(-ik 1 )Pd(ki)%(fc 1 )™ G (i? s fci)e-* 1 ' Xp(q) . 

ki 

(17) 

It is important to note that ^(q) cannot be directly 
derived with Fourier transformation because x p are not 
spaced regularly. To bypass this problem, we intro- 
duce a transitional field in Fourier space, r(ki) = 
N^(— iki)pd(k-i)Rs(ki)wG(R B ki), which is related to the 
density- vector through ^(q) = T[— x p (q)] = T[L — x p (q)], 
with r(x) the Fourier transform of r(ki) and the vector 
L = (L, L, L). One might think that the density- vector field 
can be derived straightforwardly via interpolation. Unfortu- 
nately, interpolation can cause smoothing and serious errors 
in the final estimation of the Hamiltonian force. In order to 
correct these effects, we proceed as follows. We first divide 
r(ki) by wcic(ki) to deconvolve the CIC interpolation that 
is applied later. We then Fourier transform the deconvolved 
r(ki) into real space to obtain T'(x). Finally we interpo- 
late r'(x) to the position, L — x p (q), to obtain ^(q) via 
a CIC scheme. We emphasize again that deconvolving the 
CIC kernel is crucial for obtaining an accurate estimate of 
the Hamiltonian force (see Section [3}. 

With Vl'(q) obtained, we can rewrite the Hamiltonian 
force as: 



^(k) 



E 

k 2 

E 



R z (k 2 )d5{k 2 ) 

k\ a*j(k) 

R4k 2 ) d6(k 2 ) 
'^(k) 



*MEiv^ k2q *(<o 



(ika-)**(k a ) 



(18) 



where \J , (k 2 ) is the Fourier transform of 'l'(q). Again, the 
last equation cannot be obtained directly from Fourier trans- 
forms since the signs in the exponents are not negative; it 
is obtained using the fact that e ,:k2 ' q = e " 4(_k2) q and that 
¥(-ka) = *-(k a ). 

Since d8{k 2 ) / d5j (k) is nonzero only when k 2 = ±k, the 



likelihood term of the Hamiltonian force for the real part of 
5(k) can be obtained as 



F (k) = ^) k .* 1 (k) ) 



k 2 

and for the imaginary part as 

2R z (k) 



Fr(k) 



k 2 



k-* (k) : 



(19) 



(20) 



where ^o(k) and vE'i(k) are the real and imaginary parts of 
vE'(k), respectively. It is interesting to note that the formu- 
lae for the Hamiltonian force are very similar to the gravita- 
tional force equation in Fourier space if we consider ^*(k) 
as the mass density field. 

2.4 Hamiltonian Mass and other Adjustable 
Parameters 

The method described above has two free parameters; the 
Hamiltonian masses 77lj(k) and the pseudo time T. The ef- 
ficiency of the HMC method is strongly dependent on the 
choices for these parameters. Taylor et al. (2008) suggested 
that the Hamiltonian mass for a variable is taken to be in- 
versely proportional to the width of the target probability 
distribution but the suggestion made in Hanson (2001) was 
almost the opposite. In this paper we take a different and 
much simpler approach. 

When the Hamilton equations are evolved to a pseudo 
time T, to first order approximation the change in 8j(k) is 



A^-(k) 



(k) 



K(k)- 



25 3 (k) 



mi(k) 



(21) 



.iW*) 

To ensure the convergence of the Hamiltonian system, 
A5j(k) cannot be much larger than Sj(k). We thus re- 
quire that both the first and second terms of A5j(k) be 
of the same order as or less than J Pn n (k)/2, the root mean 
square (RMS) of Sj(k). Let us first consider the second 
term. Supposing T ~ 1, one can deduce that the mass is 
of the same order as (or less than) 8 3 ;(k)/(P Un (fc)/2) 3/2 + 
Fj (k) / y/ Pun(fe) /2. We therefore define the Hamiltonian 
mass as, 



rrij(k) = m(k) 



where (■ • -)k denotes average over the phase of k. Note that 
the first and second terms in the mass equation are actually 
the RMS of 5 J (k)/(P lin (fc)/2) 3 / 2 and F^K) / yf P lia {k) /2, re- 
spectively. Since the momentum Pj(k) follows a Gaussian 
distribution with variance mj(k), i.e. pj(k) ~ y mj(k), our 
mass definition also ensures that the first term in Eq (|21|) is 
comparable to or less than \/Py m {k')/2. For consistency, we 
set n max = 13 and choose r max around 0.1 to guarantee that 
T is of order unity. 

The quantities {Ff{k))^ vary significantly before the 
HMC chain converges, and so it is not necessary to com- 
pute the masses at every step. In practice we only calculate 
the masses twice during the whole sampling. The first cal- 
culation is before the generation of the first sample. After 
proceeding iV m accepted chain steps, we use the new Hamil- 
tonian forces to update the mass variables and then retain 
the masses all the way to the end of the sampling. In the 
next section, we will show that the parameters r max and N m 
have no important impact on our final results. 
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2.5 Summary of Method 

Given the complicated, technical nature of the method de- 
scribed above, this subsection gives a step-by-step descrip- 
tion of the Hamiltonian Monte Carlo method used to re- 
construct the linear density field, given a final, evolved den- 
sity field, pt(x). This serves as a 'road-map' for anyone who 
wishes to implement this powerful method. 

(i) Pick a cosmology, which sets the analytical, linear 
power spectrum, Pn n (k). 

(ii) Randomly pick an initial guess for the modes 5(k) of 
the initial density field by specifying the corresponding real 
and imaginary parts. 

(iii) Pick a set of Hamiltonian masses, mj(k) using 
Eq. 

(iv) Randomly draw a set of momenta, pj (k) , from multi- 
dimensional, un-correlated Gaussian distributions with vari- 
ances mj(k). 

(v) Randomly pick values for the number of time steps, n, 
and the leapfrog time steps r, from uniform distributions in 
the range [0, n max ] and [0,r max ], respectively. In this paper, 
we set n max = 13 and r max = 0.1 unless otherwise specified. 

(vi) Integrate the Hamiltonian "equations of motion" us- 
ing the leapfrog technique (Eqs. [IQ]-[12]) for a total pseudo- 
time T = jit, starting from [5, (k),p., (k)]. The detailed nu- 
merical operations performed for each time step r are listed 
below. 

(vii) Accept the new 'state' [<$j(k),p<(k)], to which the 
system has evolved, with a probability given by Eq. @. 

(viii) Go back to step (iv) and repeat until the Markov 
Chain has converged and accumulated the required number 
of chain elements (i.e., properly samples the target proba- 
bility distribution ip given by Eq. [?]). 

Finally, we list the numerical operations performed in 
each leapfrog time step r: 

• Starting from the modes, 8j(k), use Eq. (|14|l to com- 
pute the displacement field, s(q), and move particles ini- 
tially located on a uniform rectangular grid of positions q 
to x p (q) = q + s(q). 

• Construct the MZA density field utilizing the clouds- 
in-cells (CIC) assignment method of Hockney & Eastwood 
(1981). 

• Fourier transform this density field using the Fast 
Fourier Transform method, and multiply the result with the 
density transfer function Rs{k) — exp[0.58<5 2 (fc/2)] and the 
Gaussian kernel WG(R s k). Divide the result with the Fourier 
transform of CIC kernel and obtain the modeled density 
field, pmod(k). 

• Use p m od(k) to compute pd(k) as described in Eq. (1131) . 
and use the method described in the paragraph below Eq. 
(|17p to compute the density- vector field *(q). 

• Fourier transform 'l'(q) to get vE'(k), and compute the 
likelihood term of the Hamiltonian forces using Eqs. (|19[) 
and (|20p for the real and imaginary parts, respectively. 

• Use the Hamiltonian forces to evolve the system accord- 
ing to Eqs. <nnji-G2jt- 

In our HMC method, the main purpose of using the 
Gaussian kernel is to suppress noises on the grid size L/N c 
that affect the efficiency of the HMC. We find that using 
R s < 2L/N C results in a quite low acceptance rate, and so 



we adopt R s > 3L/N C throughout the paper. The efficiency 
of our HMC method depends also on the value of p, in the 
sense that a smaller fi leads to a lower acceptance rate. In 
fact, the value of fj, affects the HMC in a similar way to 
the smoothing scale R s , as we will see in the next section. 
Moreover, in order to achieve a good performance, the mass 
variables and the pseudo-time T — tit should not be spec- 
ified independently but according to their combinations in 
Eq. (1211) . In this paper, we always choose T to be of order 
unity and derive the mass variables accordingly. Finally we 
note that our method is very fast. The computations shown 
below are all performed using one single processor (AMD 
Opteron 8380, 2.5 GHz). Each chain step takes about 21, 
222 and 2,080 seconds for iV c = 128, 256 and 512, respec- 
tively (see Table [TJ|. 



3 TEST WITH iV-BODY SIMULATIONS 

In this section we use the "Millennium Simulation" (MS, 
Springel et al. 2005) to test our method and tune the ad- 
justable parameters when necessary. This simulation as- 
sumes a spatially-flat ACDM model, with density parameter 
fi m = 0.25, baryon density parameter fib = 0.045, Hubble 
constant h = 0.73, and the linear RMS density fluctuation 
in a sphere of 8/i _1 Mpc radius, erg = 0.9. The CDM den- 
sity field of this simulation was traced by 2160 3 particles, 
each having a mass of 8.6 x 10 8 /i _1 Mq, in a cubic box of 
L — 500 ft _1 Mpc. We divide the simulation box into N'i grid 
cells and use a Gaussian kernel with a smoothing scale of R B 
to smooth the particle distribution on to the grid (see Ta- 
ble [T] for the values of N c and i? s used) . The method used 
to sample the density field on the grid is the same as that 
used to calculate p m od(x) except that now we do not in- 
clude the density transfer function. The resultant density 
field, denoted by Pf(x), is what we want to match in the 
reconstruction, and we assume erf(x) = ppf(x) as discussed 
in H2.1\ Moreover w(x) is always set to be unity and the 
non-linear scale fcjvL used in the transfer function for the 
Zel'dovich approximation is chosen to be 0.28 hMpc~ . 

Before showing the test results, we verify the accuracy 
of our estimation of the Hamiltonian forces. To this end, let 
us start with how the forces should be calculated based on 
their definitions. Suppose we want to calculate the Hamil- 
tonian force for a chosen variable Sj(k). We alter Sj(k) by 
a small amount, A<5j(k), with all other variables held fixed. 
This leads to a small variation, A\ 2 , in the parameter \ 2 - 
Consequently we can obtain the corresponding force numer- 
ically, Ff(k) = Ax 2 /A5 j (k). This is what we would like 
to have. Unfortunately, this method is very time-consuming 
and cannot be used to evolve the Hamiltonian system. The 
left panel of Fig. [1] shows a plot of fj(k), obtained using 
Eqs. (fT9)l and (f20)l . versus ^"(k). There is almost no visible 
difference between these two quantities. We also present a 
probability distribution of Fj :(k)/.Fj l (k) in the right panel of 
Fig. [1] The distribution exhibits a high peak at unity with 
a shallow but broad wing. Our further check finds that the 
Hamiltonian forces in the broad wing are very small, which 
explains why we cannot see any scatter in the left panel. We 
also show the Hamiltonian forces without deconvolving the 
CIC kernel in deriving the density- vector field ^(q). The 
resultant Fj(k) is systematically smaller than F™(k) due to 
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the smoothing with the CIC kernel. This demonstrates that 
our estimates of the Hamiltonian forces based on Eqs. (|19p 
and (|20p are accurate, and that deconvolution of the CIC 
kernel is essential. 

Now we apply our Hamiltonian method to the MS sim- 
ulation. We first perform a test with the following parame- 
ters: TVc = 128, i? s = 11.7/i _1 Mpc, ,u = 0.5 and r max = 0.1. 
In what follows we refer to this as our "primary test" . The 
initial set of Sj (k) is randomly drawn from the prior proba- 
bility distribution given by Eq.(fTJ. We note that our results 
shown in this paper are not sensitive to the choice of the ini- 
tial set of 5j(k). We make calculations of 2,000 HMC steps 
and the average acceptance rate is A r = 72%. Fig[2] shows 
Xw = xVEx^W (here, £ x w(x) = iV c 3 ) as a function of 
chain step. One can see that Xw drops sharply at the begin- 
ning (the burn- in phase), then remains almost constant after 
about 150 chain steps (the convergence phase). The density 
field of a converged chain step matches well the original in- 
put density field, with a RMS difference between the two 
about fj,-\/2Xw — 4.6%. For reference, the important param- 
eters and characteristics of the primary test are listed in 
Tabled 

In Fig. [3] we show the power spectra measured from 
different chain steps. We refer these power spectra as the 
Hamiltonian power spectra [Pu(k)], to distinguish them 
from the input linear power spectrum and the analytic lin- 
ear power spectrum. The 'two-phase' behavior can be clearly 
seen in the power spectra evolution. For the first 150 steps, 
one sees an obvious tuning process: on large scales the 
Hamiltonian spectrum first decreases and then increases, 
while on small scales the behavior is the opposite. The first 
(starting) spectrum is very similar to the analytic linear 
power spectrum, because the initial 5j(k) are drawn from 
Eq.{T}. After about 150 steps, the Hamiltonian spectra set- 
tle close to the initial spectrum used in the simulation, im- 
plying the convergence of the chain. Both the x 2 an d the 
power spectrum results demonstrate that our method en- 
sures quick convergence. 

Such two-phase behavior is common in Hamiltonian 
sampling (e.g. Hanson 2001; Taylor et al. 2008) and can 
be understood in terms of the behavior of a physical sys- 
tem. The kinetic energy of the system is directly related to 
the momenta, and is reset before each chain step, so that 
the kinetic energy does not change dramatically during the 
random walk. Meanwhile the system initially has very large 
potential, as well as large x 2 , as its initial position, i.e. Sj(k), 
is randomly generated. At the beginning, the system is un- 
stable and tends to fall into the deep potential well rapidly, 
which leads to decreases in x 2 an d in the potential. When the 
balance between the kinetic energy and potential is achieved 
(i.e. the burn-in phase is complete), the system becomes 
quasi-static and so x 2 remains more or less constant. In the 
convergence phase, the accepted states should all obey the 
target probability distribution, which is what we want. 

Recall again that we want to generate a linear den- 
sity field, which obeys the prior Gaussian probability func- 
tion specified by a linear power spectrum and evolves to a 
non-linear density field that matches the input density field. 
In what follows, we will examine our results in two differ- 
ent aspects. One is whether or not our reconstructed <5,(k) 
matches the prior constraints. The other is how well the pre- 



dicted (or modeled) density field (from the MZA), p mo d, and 
the original simulated density field, pf , match each other. 

The very small Xw f°r the converged chain states 
demonstrates that our method can recover the input density 
field at high accuracy. To further quantify at which scales 
Pmod matches pf well, we measure the phase correlation be- 
tween the two density fields. We define a phase correlation 
coefficient between two fields X(k) and Y(k) as 



C p (k) = X ® Y - 



(x (k)y (k) + x 1 (k)y 1 (k)) k 



(23) 



where the subscripts and 1 indicate the real and imagi- 
nary parts, respectively. Note that C p (k) — 1 means that 
the two quantities have exactly the same phase, while 
C p (k) = indicates no correlation. We show the phase 
correlation between pi and p mo( j for our primary test in 
Fig. [4] The correlation coefficient is close to one at large 
scales (k < 0.17 ftMpc -1 ), and declines quickly towards 
smaller scales. Such a sharp transition is due to the fact 
that small-scale structures are severely suppressed by our 
Gaussian smoothing and so contribute little to x 2 an d to the 
Hamiltonian force. As a result, the transition scale should 
be related to the smoothing scale. To demonstrate this, we 
perform two additional tests with R B = 23.4/i _1 Mpc and 
R s — 5.86/i _1 Mpc, assuming the same JV m , r max and jj, as 
in the primary test. We adopt iV c = 64 for the test with 
large Rs, and N c = 256 for the one with smaller Rs one. For 
R B = 5.86ft _1 Mpc the correlation coefficient remains about 
unity until k > 0.3/iMpc" 1 , while for R B = 23.4/i _1 Mpc 
the transition starts already from k ~ 0.08 ftMpc . As ex- 
pected, the transition scale is inversely proportional to the 
value of R B . We also show Xw as a function of chain steps for 
these two tests in Fig. [5] Similar to the primary test, their 
Xw first declines rapidly and then stabilizes at some small 
values. 

Now let us check whether our method can recover 
the power spectrum and the Gaussian distribution used 
as priors in our reconstruction. Inspecting the Hamilto- 
nian spectra in more detail, we see that, over the en- 
tire range of wavenumbers, the Hamiltonian spectra con- 
verge to the linear power spectrum used in the MS sim- 
ulation. However, there is noticeable discrepancy between 
the Hamiltonian spectra and the linear power spectrum on 
intermediate scales (~ 0.2/iMpc -1 ). The distributions of 
Sj (k) / Pun (k) /2 are shown in Fig. [5] for three different 
(large, intermediate and small) scales. These distributions 
can all be well fitted with a Gaussian, as expected, with 
the values of a shown in the corresponding panel. The best- 
fitting a on both small and large scales are very close to the 
expectation value, a = 1, indicating the spectrum is well re- 
produced on these scales. On intermediate scales, however, 
a deviation of a from unity is clearly seen, consistent with 
the power spectrum result. 

To understand the origin of this discrepancy on inter- 
mediate scales, we perform a series of tests with smoothing 
scales ranging from 3/i _1 Mpc to 23fo _1 Mpc. In each case, 
we find that the input linear power spectrum is well re- 
covered at both large and small scales, but that noticeable 
discrepancies are apparent on intermediate scales. The dis- 
crepancy moves gradually from small to large wavenumbers 
as the smoothing scale decreases (see e.g. the lower right 
panel of Fig. [6}. We define a wavenumber k c = ~^2 k k(l — 
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Table 1. The important parameters and characteristics for HMC and resimulation. Here L and R a are in unit of h~ 1 ~Mpc, while the 
softening length e is in unit of h~ kpc. The mean consumption time t c for each chain step is in second. A r is the acceptance rate. The 
particle mass m p for resimulation is in unit of 10 10 ft _1 M0. 
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Resimulation 






L 


N c 


R s 


Xw 


tc 




N p 


m p 


e 


Zi 


Primary Test 


500 


128 


11.7 


0.0042 


21 


72% 


128 3 


414 


80 


20 


LRR 


726 


512 


5.67 


0.0029 


2080 


35% 


512 3 


19.8 


28 


30 


HRR 


181.5 


256 


2.84 


0.0038 


222 


58% 


256 3 


2.47 


15 


36 



rO))/£ fc (l - r{k)) with r(k) = P H (fc)/Piin(fc) to quantify 
this scale, and find that the k c - P s relation can be well fitted 
by k c — 1.88/i?g' 94 . Note that the phase correlation between 
p m od and pi also depends on the smoothing scale. We show 
kc so defined as the dashed lines in the phase correlation 
plot (Fig. 3]), which clearly shows that k c also characterizes 
the transition scale in the phase correlation. The question is 
why the reconstructed linear power spectrum has the devia- 
tion from the analytical linear power spectrum just around 
k c . According to Eq.©, the Hamiltonian force has two com- 
ponents, the prior term and the likelihood term. The phase 
can be well recovered only when the likelihood term dom- 
inates the force because the prior term actually generates 
random phases. The mean ratio of these two terms is, 



R F (k) 



E J (^ 2 (k)> k 



(24) 



We find that Rp(k) decreases continuously with k and is 
about unity around k c . This implies that the \ 2 is more 
sensitive to <5(k) at smaller k, and is almost completely in- 
dependent of S(k) at k 3> k c . At k 3> k c , the trajectories 
of <5(k) in the HMC are dominated by the prior so that 
the Hamiltonian spectra match the linear power spectrum 
but the Fourier phases are not constrained. At k <C fe c , on 
the other hand, the trajectories of <5(k) is governed by the 
likelihood term so that they try to trace the original linear 
density field of the MS simulation, consequently leading to 
a small x 2 and a strong phase correlation between the re- 
constructed and original density fields. However, on scales 
around k c , 5(k) has a reduced contribution to x 2 and t ne 
posterior distribution is partly regulated by the prior. Con- 
sequently the final result is a compromise between the prior 
constraint and the likelihood. Since the likelihood term de- 
creases with increasing smoothing scale [see Eqs. (|17[l . (|19p 
and (|20p ] while the prior term is not, it explains why both 
fee and the transition scale of the phase correlation depend 
linearly on the smoothing scale. 

Although the discrepancy between the reconstructed 
and original linear power spectra is not big, we may want to 
correct it using a 'renormalization' process. First we visu- 
ally identify the region where the discrepancy is significant 
(0.12 < k < 0.21 /iM pc' 1 for the primary test), and then 
use \/ Pii n (fc) /Pn (fe) to scale the amplitude of 5j(k) with- 
out changing its phase. Since Sj(k) in this region are still 
well described by a Gaussian distribution (see Fig. [5]), and 
since rescaling preserves the shape of the distribution, the 
distribution of the scaled 5j(k) is still close to Gaussian. Be- 
cause the discrepancy to be corrected is fairly small so that 
the contribution of <5(k) around fc c to \ 2 is n °t large, the 



scaling does not cause any significant change in x and the 
phase correlation. In Section [5] we will use the renormalized 
<5(k) to generate initial conditions at high redshift, then use 
a A r -body simulation to evolve the initial condition to the 
present day and compare the resimulated density fields with 
the original simulated density field. 

Finally, we discuss the effects of changing r max , 7V m 
and p on our results. We perform several tests with dif- 
ferent Tmax, N m and p. Similar to the primary test, all these 
tests exhibit a two-phase behavior in the HMC, and a final 
convergence with a low Xw ls always achieved (see Fig. [2}. 
Inspecting the results in detail, one can see that the val- 
ues of r max and iV m affect the number of steps required for 
burn-in: large T max and small 7V m both result in a quick 
burn-in phase. Similar to R s , the value of p also affects the 
difference between the converged p mo d and pf, in the sense 
that a smaller p results in smaller py / 2x 2 D , as shown in the 
lower left panel of Fig. [5] Furthermore, the value of k c de- 
creases with the increase of p (see the lower left panel of 
Fig. [6)l , because a larger p suppresses more the contribution 
of small-scale structures to x 2 ■ Thus p affects the accuracy 
of the final result in a way quite similar to P a . On the other 
hand, since the MZA becomes increasingly inaccurate on 
small scales, a smaller p also leads to a lower acceptance 
rate. As a compromise between efficiency and accuracy, and 
because the effects of changes in p and Rs are degenerate, 
we fix p — 0.5 and only test the impact of changing P s . We 
have also checked the Hamiltonian power spectra and the 
probability distribution of <5j(k)/ *J Pu n (k) /2 in these tests 
and found that they are very similar to those in the primary 
test (Fig.©. 



4 APPLICATION TO RECONSTRUCTED 
DENSITY FIELD 

The HMC method discussed above needs a present-day den- 
sity field as an input to reconstruct the initial linear density 
field. Therefore the method is useful only when we have a 
reliable method to obtain the present-day density field from 
observation. In Wang et al. (2009), we developed a method 
to reconstruct the present-day density field from the distri- 
bution of galaxy groups (halos). In this section, we first use 
this same method to reconstruct the density field from a 
mock group catalogue and then apply our HMC method de- 
scribed above to the reconstructed density field to obtain the 
initial linear density field. Here we briefly describe the con- 
struction of the mock galaxy and group catalogues and the 
method to correct for redshift space distortions due to the 
peculiar velocities of the groups. The two mock catalogues 
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are exactly the same as those used in Wang et al. (2012) 
and are constructed from the MS simulation with the use 
of the SDSS DR7 sky coverage and selection functions. We 
refer the interested readers to Wang et al. (2009, 2012) and 
references therein for further details. 

4.1 The Mock Galaxy and Group Catalogues 

The construction of the mock galaxy catalogue is similar to 
that in Yang et al. (2007, hereafter Y07). First, we stack 
3x3x3 replicates of the MS simulation box and pop- 
ulate the dark matter haloes in the stacked boxes with 
galaxies of different luminosities, using the conditional lu- 
minosity function (CLF, Yang, Mo & van den Bosch 2003) 
model of van den Bosch et al. (2007). These dark matter 
halos are identified using the standard friends-of-friends al- 
gorithm (hereafter FOF, Davis et al. 1985) with a linking 
length of 0.2 times the mean inter-particle separation. The 
CLF describes the halo occupation statistics of SDSS galax- 
ies, and accurately matches the SDSS luminosity function 
and the clustering properties of SDSS galaxies as functions 
of luminosity. Phase-space parameters are assigned to the 
galaxies following More et al. (2009; see also Yang et al. 
2004). Briefly, the brightest (central) galaxy in each halo is 
located at rest at the halo center, while the other galax- 
ies (satellites) are distributed spherically following an NFW 
(Navarro, Frenk & White 1997) number density distribu- 
tion with the concentration-mass relation given by Maccio 
et al. (2007). At the assigned position of each satellite galaxy, 
one-dimensional velocities are drawn from a Gaussian with 
a dispersion computed from the Jeans equation assuming 
isotropy. Next we place a virtual observer at the center of 
the stacked boxes and assign each galaxy a redshift and 
(q.j, Sj) coordinates. The redshift is a combination of the 
galaxy's cosmological redshift and its peculiar velocity along 
the line-of-sight. Then we construct a mock galaxy catalogue 
by mimicking the sky coverage of the SDSS DR7, taking de- 
tailed account of the angular variations in the magnitude 
limits and the survey completeness (see Li et al. 2007 for 
details). DR7 covers two sky regions: a larger region in the 
Northern Galactic Cap (NGC) and a smaller region in the 
Southern Galactic Cap. In this paper we only use the more 
contiguous NGC region. 

Mock galaxy groups are selected using the adaptive 
halo-based group finder developed by Yang et al. (2005). 
The application of this group finder to the SDSS DR4 is 
described in detail in Y07. The application to our mock 
galaxy catalogue is exactly the same, except that the sky 
coverage is significantly larger and the adopted cosmology 
is different. The geometry of the SDSS used for the group 
catalogue is defined as the region on the sky that satisfies 
the redshift completeness criterion (C z > 0.7). As described 
in detail in Y07, our group finder takes account of the sur- 
vey edges in the SDSS volume by estimating the fraction, 
/edge, of each group's volume that falls inside the survey 
volume. Group luminosities and masses are then corrected 
for this fraction, and groups with /edge < 0.6 are excluded, 
which removes only about 1.6% of all groups. The major- 
ity of the groups in the catalogue have two estimates of 
their dark matter halo masses: one based on the ranking 
of the total characteristic luminosities of groups, and the 
other based on the ranking of the total characteristic stellar 



masses, both determined from group member galaxies more 
luminous than M r — 5 log h = —19.5. As shown in Y07, both 
halo masses agree very well with each other, with average 
scatter that decreases from ~ 0.1 dex at the low mass end to 
~ 0.05 at the massive end. In this paper we adopt the halo 
masses based on the characteristic luminosity ranking. We 
use all groups with assigned mass Mh > M t h = 10 12 /i _1 M Q 
and we restrict our reconstruction of the present-day den- 
sity field to the nearby volume covering the redshift range 
0.01 < 2 < 0.12, which we call the survey volume. 

In order to reconstruct the density field in real space, 
we have to correct for redshift space distortions. To do that, 
we follow exactly the same procedure as described in Wang 
et al. (2012, see also Wang et al. 2009). We first embed the 
survey volume in a periodic, cubic box of 726 h~ Mpc on a 
side. In the following, this will be referred to as the survey 
box to distinguish it from the survey volume and the sim- 
ulation box. We divide the box into 1024 3 grids, and sort 
them into two types: survey grids and non-survey grids, de- 
pending on whether or not the center of the grid in question 
is located inside the survey volume. We assign the mass of 
each group (with mass > M t h) on the survey grids accord- 
ing to its redshift-space coordinates. Non-survey grids are 
assigned a density equal to the average mass density of the 
groups (with Mh > M t h) in the survey volume so that the 
entire survey box has the same mean density as the actual 
survey volume. Then we compute the overdensity field of our 
groups in redshift space, and smooth it using a Gaussian ker- 
nel with a smoothing scale of 8 h~ 1 Mpc. In the linear regime, 
the peculiar velocities induced by density perturbations are 
proportional to the amplitudes of the density fluctuations, 
hence we can use this smoothed overdensity field to derive 
the peculiar velocities. Adopting the relatively large smooth- 
ing scale can effectively suppress non-linear velocities that 
cannot be predicted accurately with the linear model. As 
shown in Wang et al. (2009), the resultant velocities based 
on the halo (group) population are biased but tightly pro- 
portional to the real velocities. We thus can predict the pe- 
culiar velocity of each group by simply taking into account 
the bias factor of the density field represented by the groups, 
and we use equation (10) in Wang et al. 2009 to calculate 
the bias factor. Finally we use the line-of-sight component 
of the predicted velocity to correct the redshift space distor- 
tions. Since the velocity field is computed using the redshift 
space distribution of the groups, this method needs to be 
iterated (typically twice) to make sure that convergence is 
achieved (Wang et al. 2012). 

4.2 Density Reconstruction from Mock Group 
Catalogue 

In Wang et al. (2009), we developed a method to reconstruct 
the present-day cosmic density field starting from the dis- 
tribution of galaxy groups (haloes). We tested this method 
using FOF haloes but did not apply it to mock groups to 
examine its reliability against uncertainties arising from sur- 
vey boundaries and false identifications of group members 
by our group finder. Survey boundaries have no direct im- 
pact on our reconstruction, but can affect the correction for 
redshift space distortions. In this subsection, we reconstruct 
the cosmic density field from the mock group catalogue cor- 
rected for redshift space distortions as described above and 
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compare our reconstructed density field with the original 
simulated density field. 

As in Wang et al. (2009), we first calculate the density 
profiles in 'domains' using the MS simulation. For a given 
population of FOF haloes with mass of Mh > Mth, we par- 
tition space into a set of domains, each of which contains 
one halo. The domain of any halo is defined in such a way 
that each point in the domain is closer (based on a distance 
measure to be defined below) to the halo than to any other 
haloes in the halo population. The proximity of a point to 
a halo with viral radius of Rh is defined as, 

r P = , (25) 

where rn is the physical distance between the halo center and 
the point. We calculate the cross correlation functions be- 
tween haloes and dark matter particles in their correspond- 
ing domains, i.e. the average density profile within the do- 
main. Fig. [7] shows these density profiles for haloes in six 
mass bins, using M th = 10 12 fc _1 M Q . Despite using a differ- 
ent cosmology, these profiles are very similar to those shown 
in Wang et al. (2009). Then we 'convolve' our mock groups 
with these density profiles to reconstruct the cosmic density 
in the following way. For a mock group of mass above the 
mass threshold Mth, we pick a density profile shown in Fig. 
\7\ according to the mass of this group. We use a Monte Carlo 
method to put particles around this group up to ~ 30 times 
the virial radius (sufficient to cover the domain of the group) 
regardless of the domain. We remove particles outside of the 
domain of the group. We repeat the above three steps for all 
groups with mass larger than M t h- Eventually we obtain a 
reconstructed density field sampled by particles, embedded 
in the survey box. 

One advantage of our reconstruction method is that we 
define a very special cross correlation function. It is different 
from the conventional cross correlation function, which is de- 
fined by the mass density within spherical shells of a given 
radius centered on a halo, regardless whether or not the 
particles are in the domain of the halo. In the conventional 
one, all haloes contribute to the average density profiles at 
all scale, especially on large scales. This has the effect of 
smoothing halo masses over very large scales. On the other 
hand, in our cross correlation based on domains, massive 
haloes (Mh > M t h) contribute only to the density profiles 
within the virial radii, while the density profiles at larger 
scales are contributed by diffuse mass and low-mass haloes 
with Mh < M t h- Another advantage of our method is that 
we use galaxy groups (haloes) instead of galaxies. While the 
bias of galaxy distribution relative to the underlying den- 
sity field may depend on various galaxy properties, such as 
luminosity and color, and the exact form of the bias is not 
well established, the use of galaxy groups (haloes) automat- 
ically takes into account the bias of the galaxy distribution 
through the connection between galaxies and dark matter 
haloes. Moreover, the usage of groups can effectively miti- 
gate the problem due to the finger-of-god effect. 

In order to compare our reconstructed density field with 
the simulated density field, we divide the survey box into 
1024 3 grids with size of 0.71 h Mpc, and smooth the sam- 
pled particles on the grids using a Gaussian kernel with a 
smoothing radius of 2.84 h~ Mpc. Since the peculiar veloci- 
ties are predicted more accurately in the inner region of the 
survey volume than near the survey boundary (Wang et al. 



2012) and we use these peculiar velocities to correct for the 
redshift space distortions, the reconstruction is expected to 
be more reliable in the inner region. To check whether or 
not this is the case, we need to define the distance of a grid 
to the boundary. Following Wang et al. (2012), we intro- 
duce a parameter, the filling factor p/, to characterize the 
closeness of a survey grid to the boundary. For a grid, g, the 
filling factor is defined as the fraction of grids in a spheri- 
cal volume of radius 80 h~ Mpc centered on the grid g, that 
are located within the survey volume. Hence, pt is expected 
to be much less than unity for a grid located close to the 
boundary, while pt ~ 1 for a grid that is located more than 
80/i _1 Mpc from any survey boundary. 

In Fig.|8]we compare the simulated density field pt and 
the reconstructed density field p IC in the inner region with 
pt > 0.6, which is about 66% of the survey volume. The 
solid line represents the mean relation and the error bars 
indicate the standard deviation in pt for a given p Ic . The 
bias in the mean relation is very small and the uncertainties 
are about 30% to 50% of p TC in most bins. At the two largest 
density bins, there is a significant deviation from the one to 
one relation. The volume of the grids in these two bins is 
tiny and their presence does not affect our reconstruction 
significantly. We then show the same comparison for grids 
near the boundary, i.e. with pt < 0.6. As one can see, the 
result is as good as that in the inner region. Apparently, 
the smoothing used is able to remedy partly the problem 
in the correction of the redshift space distortions near the 
survey boundary. Overall, with an appropriate choice of the 
smoothing radius, our method is able to reconstruct the den- 
sity field accurately, and the effects due to survey boundary 
and group contamination do not introduce any significant 
bias in the reconstruction. 

In the next subsection, we will apply our HMC method 
to the reconstructed density field in two volumes, one is the 
entire survey volume and the other is a cubic volume well 
inside the survey volume (see below for details). Before do- 
ing that, we show in Fig[5] the phase correlation between 
our reconstructed and simulation fields in these two vol- 
umes. These correlations set an upper limit on the accuracy 
we can achieve for our reconstruction of the linear density 
field. As one can see, the correlation coefficient is almost 
unity on large scales, and declines gradually with the in- 
crease of the wavenumber. At k ~ 1/iMpc -1 the coefficient 
drops to ~ 0.5, demonstrating that our method is success- 
ful well into the nonlinear regime (the translinear scale is 
k ~ 0.15 /iMpc -1 ). The phase correlation at k > 0.2/iMpc" 1 
obtained for the inner cube is slightly stronger than that for 
the total volume, indicating that boundary effects indeed 
have a non-negligible impact on the reconstruction. 

4.3 Reconstruction of Linear Density Field 

We first apply our HMC method to the reconstructed den- 
sity field in the entire survey volume, which is embedded in 
the survey box with size L = 726 ft _1 Mpc on a side. In order 
to examine the ability of our method over a large dynamical 
range, we divide the survey box into iV c = 512 grids in each 
dimension, so that our reconstruction deals with more than 
Nc ~ 10 s free parameters. We adopt a smoothing radius 
of Rs = 5.67fe _1 Mpc, on which our reconstructed density 
field is quite similar to the original density field. In order 
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to derive the weight w(x), we divide each grid into 2 3 sub- 
grids. If more than 6 subgrids of a grid are located inside the 
survey volume, this grid is assigned a weight of unity, oth- 
erwise zero. The other parameters are chosen as r max =0.1, 
/i = 0.5 and iV m = 50, similar to those in the primary test 
described above. This application is referred to as the 'Low 
Resolution Run' or LRR in the following. As shown in Ta- 
ble [1] on average it takes ~ 2080 seconds for each HMC 
step and the acceptance rate is about 35%. The Xw value 
for a converged state is about 0.0029, indicating that the 
RMS difference between p mo d(x) and p rc (x) is only 3.8%. 
We show the phase correlation between p mo d and p rc in the 
left panel of Fig. [9] The correlation is almost unity at large 
scales, and drops quickly around k c , consistent with our test 
results presented above. 

In the left panel of Fig. [10] we compare a converged 
Hamiltonian spectrum, Pu(k), with the analytic linear 
power spectrum, Pn n (k), used in the prior. As one can 
see, Pn(k) is very close to Pii n (fc) over the entire range 
of scale. In particular, the discrepancy around k c is also 
small. The distribution of the reconstructed initial density 
field is extremely close to Gaussian, as shown in the up- 
per three panels of Fig. 1111 As described above, we can 
also renormalize the reconstructed linear power spectrum 
by visually identifying the discrepancy region (0.28 < k < 
0.47 ftMpc -1 ) and scaling the corresponding Sj(k) with a 
factor v/ Pijn (k) I Ph (fe) . The new power spectrum so ob- 
tained and the distributions of the re-normalized <5j(k) are 
also shown in the corresponding figures. Since the discrep- 
ancy is tiny, the curves before and after the renormaliza- 
tion are undistinguishable in the figure. As shown in the 
upper panels of Fig. 111! these re-normalized 5j (k) are well 
described by a Gaussian function with unity variance. These 
results demonstrate clearly that our HMC method works 
well on large scales. 

To examine the performance of our HMC method on 
small scales, we perform a High Resolution Run (HRR) with 
a small smoothing scale. This is a cubic box of 100 h~ Mpc 
on a side located in the inner region of the survey volume, 
put inside a larger periodic box with L = 181.5 fe _1 Mpc. We 
divide the larger box into N c — 256 grids in each dimension, 
and adopt a smoothing scale of R B = 2.84 h~ 1 Mpc, for which 
our reconstructed density field is in good agreement with the 
original simulation. Only grid cells that are located within 
both the small box and the survey volume are assigned a 
weight of unity. All other grid cells are assigne a weight of 
zero. Note that some grids in the small box may not be in 
the survey volume, because of the existence of small holes in 
the survey mask. The other parameters are chosen to be the 
same as in the LRR. On average the HRR takes 222 seconds 
for each HMC step, and the acceptance rate is about 58%. 
The chain finally converges to Xw — 0.0038, corresponding 
to a RMS of 4.4% in the difference between p mo d(x) and 
p rc (x). The corresponding phase correlation is shown in the 
right panel of Fig. [5] Here we see a significant correlation 
all the way to k ~ 1 ZiMpc^ 1 , a scale much smaller than the 
translinear scale which corresponds to k ~ 0.15 /iMpc -1 . 

The converged Ph(&) and the distributions of 
(5j(k)/ 1 /Pi in (fc)/2 for the HRR are shown in Figs. [TOl and 
llll Here we see a significant bump at 0.3 < k < 1 /i _1 Mpc 
Pn(k) compared with Pn n (k). The reason for this bump is 
twofold. First, the bump is around k c , suggesting that it 



is a compromise between the prior constraint and the like- 
lihood (see Section [3]). Second, the bump is partly due to 
the inaccuracy in the adopted approximation of the struc- 
ture formation model on small scales. According to TZ12, 
the mode-coupling term, pA/c(k), cannot be neglected at 
k > 0.4/iMpc -1 (see their figure 1). On such scales, the 
amplitude of p mo d(k) predicted by the MZA is somewhat 
lower than that of the fully evolved density field. To achieve 
a small x 2 m the HMC, the Hamiltonian spectra have to 
be enhanced to compensate the deficit. On even smaller 
scales, k > lftMpc -1 , however, the Hamiltonian force is 
dominated by the prior term so that Pn(k) is forced back 
to Pnn(fc). Despite of this inaccuracy, our HMC method can 
still recover more than half of the phase correlation all the 
way to k ~ 1/iMpc -1 , as demonstrated in the next sec- 
tion. As before, we re-normalized the amplitudes of the 
Fourier modes to correct for the discrepancy in the range 
of 0.35 < k < lfaMpc -1 . The corrected power spectrum 
and distribution of the re- normalized <5j(k), are shown in 
the right panel of Fig. [10] and the lower panels of Fig llll 
respectively. 



5 RE-SIMULATIONS OF RECONSTRUCTED 
LINEAR DENSITY FIELD 

Up to now, all our results are based on the structure for- 
mation model of TZ12. The advantage of using the TZ12 
model is that it is very fast and can thus be implemented 
into our HMC to infer the initial linear density field. How- 
ever, the TZ12 model, which is based on the Zel'dovich 
approximation, is not expected to work accurately in the 
highly non-linear regime. Although the modeled density field 
shown above is closely correlated with the input present-day 
density, it is unclear to which extent the density field fully 
evolved from our reconstructed linear density field matches 
the original simulated density field, especially on small scales 
where non-linear evolution becomes important. In order to 
test the model in the highly nonlinear regime, we need to use 
full JV-body simulations to evolve our reconstructed initial 
linear density field to the present day and to compare the 
re-simulated density field with the original simulated density 
field. 

To do this we set up the initial condition for our con- 
strained N-body simulation using the renormalized S(k). We 
generate the initial displacement field at a given high red- 
shift, Zi, using the Zel'dovich approximation. The displace- 
ment field is used to perturb the positions of the N-body 
particles that initially have uniform distribution. Each par- 
ticle is assigned a velocity according to the growing mode 
solution of the linear perturbations. We then use the N-body 
code Gadget-2 (Springel 2005) to evolve the initial condition 
to the present day. The fully evolved density field is referred 
to as the resimulated density field, and denoted by p rs . 

5.1 Initial Conditions from Simulated Density 
Field 

Let us first consider S(k) reconstructed from the origi- 
nal simulated density field. For our primary test, we use 
N p = 128 3 particles in a box of L — 500/i _1 Mpc to trace 
the evolution of the density field. As shown in Table. [1] we 
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adopt an initial redshift = 20, a particle mass m p ~ 4.14x 
10 12 ft _1 M0 and a force softening length e = 80 ft _1 kpc. The 
dashed black line in Fig. [4] shows the phase correlation be- 
tween p rs and p\. As can be seen, this phase correlation is 
very similar to that between p m od and pf. We also show 
the results of the resimulation from the reconstructed linear 
density field using R B — 5.86 /i -1 Mpc and 3.91 /i _1 Mpc, re- 
spectively. For the test with R B = 5.86 h~ 1 Mpc, the phase 
correlation is again very similar to that between p mo d and pf , 
suggesting that the MZA works well up to k Si 0.3 /iMpc -1 . 
In the test with R s = 3.91 h~ Mpc, however, the p rs - p{ 
correlation is much stronger than the p mc ,d - Pf correlation 
Sit k kc . 

Given that the HMC method tends to minimize the dif- 
ference between p mo d and p rc , it is unexpected that the p mo d 
- pf correlation is much weaker than the p rs - pi correlation 
at k > k c . The most important difference between p mo d and 
p rs is that the later is the result of fully non-linear process, 
in which the mode of the non-linear density field on small 
scales may be coupled to that on large scales (see, e.g., Tas- 
sev & Zaldarriaga 2011). Consequently, part of phases at 
k > k c , where the initial phases are not well constrained, 
is reproduced in the re-simulation. Such mode coupling is 
not included in the TZ12 model based on quasi-linear the- 
ory, so that the phase information in pf at k 3> k c is almost 
completely lost in p mo d- 

Our above results clearly show that the match between 
the resimulation and the original simulation is better at 
smaller R s . To quantify this trend, we introduce a quan- 
tity fch to measure the scale at which p rs can well match 
Pf. It is defined in such a way that the phase correlation 
between the two density fields at kh is half. In Fig. 1121 we 
show fch as a function of R s . On large scales, kh increases 
with decreasing Rs, consistent with expectation because in- 
formation of the density field on scales below Rs is lost due 
to the smoothing. As it reaches about unity, however, the 
value of kh becomes insensitive to R s . This is also expected, 
because the TZ12 model is not expected to work accurately 
on very small scales. This demonstrates that in order to re- 
construct the density field on scales below k ~ 1/iMpc -1 , 
one has to use a model that is more accurate than the TZ12 
model adopted here. 

5.2 Initial Conditions from Reconstructed 
Density Field 

In this subsection, we use JV-body simulations to evolve the 
initial conditions obtained from the reconstructed density 
field (Section [4j. For reference we list the parameters for 
both LRR and HRR in Table [T] To inspect our results visu- 
ally, we present density maps of the same thin slices in the 
original simulation used to construct the mock catalog, the 
density field reconstructed from the mock catalog, and the 
density field in the re-simulations. Fig. H3l shows the results 
for the LRR. Here all density fields are smoothed with a 
Gaussian kernel with R s — 5.67 /i -1 Mpc. Within the survey 
volume, the three maps look quite similar; almost all struc- 
tures in the original simulation, such as massive clusters, 
filaments and underdense voids, are reproduced in the re- 
simulation. The density maps for the HRR, smoothed with 
R s = 2.84/i _1 Mpc, are presented in Fig. 1141 There are about 
twenty bright spots in the original simulation, which corre- 



spond to a single halo or a cluster of a few haloes with 
masses down to about a few times 10 12 /i _1 M Q . Most of 
these structures are clearly reproduced in our resimulation. 
In addition, some small filaments in the original simulation 
are also correctly reproduced in the re-simulation. This is 
quite remarkable, given that our reconstruction from the 
mock catalog uses only groups with assigned masses above 
Mth = 10 12 Ii~ 1 Mq and that non-linear effects are impor- 
tant on such small scales. 

In the left panel of Fig. [15] we show the comparison of 
p rs with p rc and pt for the LRR. There is weak bias at the 
high density end in both relations: while p rs is higher than 
p rc , it is slightly lower than pt . Since the initial condition for 
the resimulation is constrained by the reconstructed density 
field, the p rs - Pre relation is much tighter than the p rs - 
pt relation. The typical dispersion in the former relation is 
about 0.05 dex, while it is about two to three times larger in 
the latter relation. The phase correlations among the three 
density fields are shown in the left panel of Fig. [9] Similar 
to the p mo d - Pre correlation, there is a sharp transition in 
the correlation coefficient from unity to about zero around 
k c . Upon closer examination, we find that the p rs - pi phase 
correlation is always lower than the minimum of the p rs - 
pr C and pre - pt phase correlations. This is expected, as the 
accuracy of the resimulation depends both on the accuracy 
of the HMC method, which determines the strength of the 
Prs-Prc correlation, and the accuracy in the reconstruction of 
the present-day density field, which determines the strength 
of the prc-pi correlation. 

The comparisons of p rs with p rc and pt for the HRR are 
shown in the right panel of Fig. [TSJ As one can see, p rs is 
also strongly correlated with both p rc and pf. The scatter 
in the p rs -pf correlation is less than 0.2 dex in most density 
bins. At the high-density end, a weak bias is present in the 
Prs- Pre relation, but such a bias is absent in the p rs - pf 
relation. At the low-density end (p < 0.3p), p rs is correlated 
neither with p rc nor pt . The reason is that in our reconstruc- 
tion of the present-day density field using mass distributions 
around halos, the minimum of the density profiles in the do- 
main is about 0.25p (Fig. 0. Therefore the information in 
the most underdense regions is totally lost when reconstruct- 
ing the present-day density field. This bias can be mitigated 
by adopting a smaller mass threshold Mth f° r the group cat- 
alogue. The usage of a smaller M t h can lower the minimum 
density that our reconstruction can reach (see Wang et al. 
2009). 

The phase correlations for the HRR are presented in 
the right panel of Fig. The phase correlation between 
p rs and pre declines gradually around k c . Similar to the re- 
sults shown in the previous subsection, this correlation is 
much stronger than the p m od - pic correlation at k > k c . 
As discussed above, this is due to mode coupling which is 
not fully included in the TZ12 model. Moreover, the p rc - 
pt phase correlation lies below the p rs - Pre phase correla- 
tion at almost all scales. It indicates that, on the smoothing 
scale R B = 2.84/i _1 Mpc, the accuracy of the re-simulation 
in matching the original density field is mainly limited by 
the reconstruction of the present day density rather than by 
the HMC method. Despite all of these, the match between 
the resimulation and the original simulation is remarkable. 
At k = fc c , the phase correlation between p rs and pt is still 
as high as 0.6; even at k = 1.0 ZiMpc -1 , about half (47%) of 
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the phase information is reproduced. We recall again that 
the translinear scale corresponds to k — 0.15 ftMpc -1 . 

We also perform tests with other values of R s and mea- 
sure fch, at which the phase correlation between the res- 
imulation and original simulation is half. We show kh as a 
function of R s as the dashed line in Fig. 1121 One see that the 
value of feh first increases with decreasing R s then remains 
almost constant, and has the maximum of 0.94 ftMpc -1 at 
R s = 2.84 /i _1 Mpc. The curve generally lies below that 
based on original simulated density field (solid line), be- 
cause the reconstruction of the present-day density field is 
not perfect. Note that at R s — 2.84 /i _1 Mpc, the value of fch 
obtained from the reconstructed present-day density field 
is similar to that obtained from the simulated present-day 
density field, because on such small scales the inaccuracy of 
the TZ12 model becomes the dominating source of error in 
the reconstructed initial condition. 



6 SUMMARY AND DISCUSSION 

In this paper we have developed a new, effective method to 
reconstruct the linear density field that underlies the for- 
mation of the cosmic density field in the local Universe. To 
this end we have developed a Hamiltonian Markov Chain 
Monte Carlo (HMC) method which allows us to sample the 
linear density field based on a posterior probability function. 
This distribution consists of two components, a prior term 
which takes into account the Gaussianity and power spec- 
trum assumed for the linear density field, and a likelihood 
term designed to ensure that the predicted density field from 
the initial condition matches a given final density field. We 
adopt the modified Zel'dovich approximation developed by 
TZ12 to model the final, evolved density field off the initial, 
linear density field. 

The HMC method is based on an analogy to a phys- 
ical system. The potential is taken to be minus the loga- 
rithm of the posterior function. The momenta are drawn 
from given Gaussian distributions before each chain step so 
that the fictitious system can continue to equilibrate and 
"orbit" within the extended potential well with the pas- 
sage of "time". The system eventually converges to a state 
in which balance between kinetic and potential "energy" is 
achieved. Using a simulated density field, we demonstrate 
that our HMC method converges very quickly, and that the 
converged linear density fields closely follow a Gaussian dis- 
tribution with a spectrum that accurately matches the input 
linear power spectrum. A small discrepancy is found at the 
scales where the likelihood and prior terms in the Hamilto- 
nian force are comparable. This discrepancy, however, can 
straightforwardly be corrected for by re-normalizing the am- 
plitudes of the corresponding Fourier modes (while keeping 
their phases fixed) with the input linear spectrum. We find 
that the modeled density field matches the input density 
field with high precision, with a RMS difference typically 
smaller than 5%. 

Since our HMC method needs the present-day den- 
sity field as a constraint in reconstructing the initial linear 
density field, we also present a method to reconstruct the 
present-day density field from mock galaxy and group cat- 
alogues. The mock catalogues are constructed from the MS 
simulation for the SDSS DR7, taking detailed account of 



the angular variations in the magnitude limits and the sur- 
vey completeness (see Wang et al. 2012 and the references 
therein) so that we can verify the reliability of our method 
in real applications. We use the method developed by Wang 
et al. (2009) to reconstruct the density field based on the 
mock group catalogue, taking into account inaccuracies in 
the group finder, as well as uncertainties arising from red- 
shift space distortions and the assignment of a halo mass to 
each individual group. We find that the phase correlation 
between the reconstructed and simulated density fields is 
almost perfect at large scales, with a correlation coefficient 
close to one. The scale at which the correlation coefficient 
drops to 0.5 is k ~ 1/iMpc -1 , indicating that our method 
works surprisingly well down to scales that are well into the 
non-linear regime. 

We apply the HMC method to the reconstructed den- 
sity field in two different volumes. The first one (LRR) is 
the entire survey volume of the SDSS embedded in a pe- 
riodic box of size 726/i _1 Mpc, and the second one (HRR) 
is a cubic box of size 100/i _1 Mpc covering the inner region 
of the survey volume. These two applications are used to 
test and verify the performance of the HMC method over 
a wide dynamic range. As an additional test of the perfor- 
mance of our methods, we use the reconstructed linear den- 
sity fields of LRR and HRR to generate initial conditions, 
which we subsequently evolve to the present day using a 
TV-body simulation code. Both visual inspection and quan- 
titative analysis show that the density field obtained from 
these re-simulations accurately match the density field of the 
original simulation used to construct the mock catalog. In 
particular, the phase correlation between re-simulation and 
original simulation has a coefficient close to unity on large 
scales and only starts to drop to 0.5 at k ~ 1 hMpc~ . This 
clearly demonstrates that our HMC method together with 
the reconstruction method of Wang et al. (2009) provides a 
robust way to reconstruct the initial conditions for the local 
cosmological density field from observational data. 

Numerous studies in the past have tried to infer the ini- 
tial conditions of structure formation in the local Universe 
using observational data such as galaxy distributions and/or 
peculiar velocity surveys (e.g. Nusser & Dekel 1992; Wein- 
berg 1992; Kolatt et al. 1996; Klypin et al. 2003; Dolag et 
al. 2005; Lavaux 2010). Most of these studies integrate the 
observed density field backwards in time to some initial time. 
However, these approaches suffer from complications in the 
observational data, such as spatial variations in the magni- 
tude limit and complex survey boundaries, as well as the 
amplification of noise and numerical errors by the decaying 
mode during backward integration (Nusser & Dekel 1992). 
As pointed out by JW12, these problems can be overcome by 
the HMC method. In our HMC method, the survey geome- 
try is taken into account by including a weight field into the 
likelihood. The amplification of noise and numerical errors 
is also not an issue in the HMC method since it uses for- 
ward evolution of the cosmic density field (see also Kitaura 
2012 for another forward-evolution method). Finally, it is 
worth emphasizing that some of the previous methods used 
had to Gaussianize the inferred initial density field using 
some order-preserving transformation. In the HMC method 
adopted here, however, the initial density field is Gaussian 
by the construction of the posterior. 

Our own HMC method has some unique advantages. 
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We design the likelihood using the present-day density field 
reconstructed from galaxy groups (i.e., dark matter haloes), 
rather than the galaxy distribution itself. The latter requires 
a detailed understanding of how galaxies are biased rela- 
tive to the underlying dark matter distribution (see e.g. Ki- 
taura & Enfilin 2008). As mentioned earlier, this bias be- 
tween galaxies and dark matter is far from trivial; it de- 
pends on galaxy properties, has stochastic and non-linear 
contributions, and may even be non-local. Adopting a sim- 
ple linear bias model would significantly underestimate the 
density in high density regions. Currently it is still unclear 
how to directly use the galaxy distribution to model the den- 
sity field, especially in high density region, in an unbiased 
way. Another problem with using the galaxy distribution 
is that the constraint becomes very poor in underdense re- 
gions, where only few or no galaxies can be observed in a 
uniformly selected galaxy sample. In our approach, these 
problems are largely absent. As shown in Section 14.21 the 
reconstructed density fields based on the detailed mock cat- 
alogues match the input density field very well over a large 
dynamical range. Furthermore, since our reconstruction re- 
lies on groups (haloes), we can accurately reproduce the high 
density regions within individual haloes. In underdense re- 
gions, using the density profiles in the domains of halos can 
recover the density field down to £ 0.25p. And in principle, 
we can recover even lower densities by simply using groups 
(haloes) above a lower mass threshold, although this requires 
either a deeper redshift survey, or the use of a more lim- 
ited volume. The reliability of our reconstruction is further 
demonstrated by the fact that the resimulations from the 
initial conditions constrained by the reconstruction match 
the original simulation remarkably well in the same density 
range. 

Moreover, our HMC method works in Fourier space. 
Different Fourier modes are mutually independent, and so 
are the real and imaginary parts of individual modes. This 
enables us to derive simple formulae for both the prior and 
likelihood terms in the Hamiltonian force, and makes the 
computation much faster. As shown in Table [T] it takes, on 
average, only about 21, 220 and 2100 seconds for each step 
for iV c = 128, 256 and 512 respectively. Particularly, our 
method can successfully handle more than one hundred mil- 
lion free parameters! Such efficiency is crucial when aiming 
to achieve high resolution in a large volume. 

In forthcoming papers, we will apply our reconstruction 
and HMC methods to the SDSS DR7 group catalogue in or- 
der to generate the initial conditions for structure formation 
in the SDSS volume. We will then use these initial condi- 
tions to run constrained simulations to study the evolution 
of the local cosmic density field. This will provide a unique 
opportunity to further our understanding of the formation 
and evolution of the galaxies we directly observe. For ex- 
ample, one can investigate the correlation between the large 
scale environments, measured from the constrained simula- 
tion, and the observational galaxy properties. Recent studies 
have found that the properties of dark matter haloes depend 
significantly on their large-scale environments, in particu- 
lar the large scale tidal fields (see e.g. Wang et al. 2011), 
and it would be interesting to see whether this is also the 
case for galaxies that reside in haloes. One can also per- 
form semi-analytical models of galaxy formation using halo 
merger trees extracted directly from the constrained sim- 



ulations. The comparison between model galaxies and real 
galaxies in the same large scale structures, such as filaments, 
sheets and clusters, will provide us an avenue to constrain 
galaxy formation in a way that is free of cosmic variance. 

Finally, the constrained simulation can also be used to 
study the dynamics and physics of the intergalactic medium 
(IGM). For instance, the hot gas and peculiar velocities ex- 
pected from the structures in the SDSS DR7 predicted by 
the constrained simulations can be used to make predic- 
tions for the Sunyaev-Zel'dovich effects (both thermal and 
kinetic), which can be compared with (forthcoming) obser- 
vations from e.g. the Atacama Cosmology Telescope (ACT; 
Swetz et al. 2011) and the PLANCK satellite. Moreover, a 
comparison of the simulated density field with quasar ab- 
sorption lines in a wide range of ionization potentials can 
provide constraints on the temperature and metallicity of 
the filaments and sheets that make up the cosmic web (Cen 
& Ostriker 1999). Such studies will provide a unique way to 
understand the nature of the absorption systems at law-z 
and the implications for the state and structure of the IGM. 
In particular, it will allow a detailed exploration of the con- 
nection and interaction between the IGM and the galaxy 
population. 
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Figure 1. Comparison between the Hamiltonian forces Fo ; i(k) and i^-^k). Fo ; i(k) is calculated using Eq. II19H and I I20I I. while Fq.^Qs.) 
is obtained using a numerical method presented in Section [3] The subscript and 1 denote the real and imaginary parts respectively. 
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Figure 2. Evolution of Xw f° r various tests. See text (Section[3j for the definition of Xw- The black lines in the four panels are the same, 
the result of primary test with parameters, N c = 128, R B = lX,7h~ Mpc, /i = 0.5 and r max = 0.1. The labels in each panel show the 
parameters different from the primary test for the corresponding test. 
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Figure 4. Solid lines: the phase correlation between p mo< j and pi for tests with different smoothing scales as indicated in the figure. 
Dashed lines: the phase correlation between p TB and pf. pf is the density of original simulation, p mo( i is the modeled density field and p TB 
is the resimulated density field. The vertical dashed lines indicate the characteristic wavenumber k c (see Section \3\ for the definition). 
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Figure 5. The black lines show the probability distribution of <5o ; i (k) / sj P\i n (k)/2 of the primary test at three different scales. For 
comparison, we also show the best-fitting Gaussian curves and their a in red. 
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Figure 6. Power spectra for various tests. In all four panels, dotted lines show the analytic linear power spectrum, Pii n (fc), the black solid 
lines show the power spectrum measured from the MS simulation and the green lines show the convergent power spectrum of primary 
test. The labels in each panel show the parameters different from the primary test for the corresponding test. 
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Figure 7. The density profiles of mass in and around the haloes in various mass bins. Here the mass threshold for the halo population 
is M t h=10 12 h~ 1 M.Q. The radius r is scaled by halo virial radius R^, and the density is scaled with p, the mean density of the universe. 
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Figure 8. The comparison of density between the original simulation and the reconstruction. The reconstruction here is obtained by 
using mock groups with M t h=10 12 fe _1 Mg and density profiles shown in Fig. [7] The two density fields are smoothed with Gaussian 
kernel with smoothing scale of 2.84/i — x Mpc. The black and red lines show the results for the grids with different filling factor, pf (see 
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Figure 9. The phase correlation among modeled density p mo d, resimulated density p ra , reconstructed density p rc and original density pf . 
The left panel shows the results for LRR. The right panel shows the results for HRR. The vertical dashed lines denote the characteristic 
wavenumber, k c . 
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Figure 10. Black lines show the analytic linear power spectrum, Pii n (fc), the red lines show one convergent Hamiltonian spectrum, Pjj(fc), 
and the blue lines show corrected power spectrum, P c (k). The left and right panels show the results for LRR and HRR, respectively. 
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Figure 11. The probability distribution of <5o ; i (k) / \J P]j n (fc)/2 at large (left), intermediate (middle) and small (right) scales for LRR 
(upper) and HRR (lower). The dashed and solid black lines show the results before and after correcting for the discrepancy, respectively. 
The dashed and solid lines are undistinguishablc at all panels except the lower middle panel. The red lines are the best-fitting Gaussian 
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Figure 12. The figure show k^, at which the phase correlation between the resimulated density field and original simulated density field 
is half, as a function of R s . The solid line show the results for the resimulations constrained directly from original simulation, while the 
dashed line show the results for resimulations constrained from the reconstructed density field. 




Figure 13. The low-resolution density maps in a slice of 600 X 300 X 4.25 h _1 Mpc. The left panel shows the resimulated density map from 
LRR. The middle and right panels show the original simulated and reconstructed density fields. All these density fields are smoothed 
with Gaussian kernel with R a = 5.67/i _1 Mpc and scaled with the mean density of the universe. 
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Figure 14. The high-resolution density maps in a slice of 100 X 100 X 6.4 h~ Mpc. The left panel shows the resimulated density field from 
HRR. The middle and right panels show the original simulated and reconstructed density fields. All these density fields are smoothed 
with Gaussian kernel with Rs = 2.84 ft - x Mpc and scaled with the mean density of the universe. 
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Figure 15. Black lines show the comparison between p rs and p rc . Red lines show the comparison between p rs and pf. All these densities 
are scaled with p, the mean density of the universe. The left panel shows the results for LRR and the densities are smoothed with 
Gaussian kernel with R s = 5.67 h~ Mpc. The right panel shows the results for HRR and the densities are smoothed with Gaussian 
kernel with R B = 2.84/i _1 Mpc. The dot lines indicate the one to one relation. 



